import math

# 定义常量
RADIANS_PER_DEGREE = math.pi / 180.0
DEGREES_PER_RADIAN = 180.0 / math.pi

# WGS84 参数
WGS84_A = 6378137.0  # 长轴
WGS84_E = 0.0818191908  # 第一偏心率

# UTM 参数
UTM_K0 = 0.9996  # 比例因子
UTM_FE = 500000.0  # 假东
UTM_FN_N = 0.0  # 假北（北半球）
UTM_FN_S = 10000000.0  # 假北（南半球）
UTM_E2 = WGS84_E ** 2  # e^2


def LLtoUTM(lat, lon):
    zone_number = int((lon + 180) / 6) + 1
    if lat >= 56.0 and lat < 64.0 and lon >= 3.0 and lon < 12.0:
        zone_number = 32
    if lat >= 72.0 and lat < 84.0:
        if lon >= 0.0 and lon < 9.0:
            zone_number = 31
        elif lon >= 9.0 and lon < 21.0:
            zone_number = 33
        elif lon >= 21.0 and lon < 33.0:
            zone_number = 35
        elif lon >= 33.0 and lon < 42.0:
            zone_number = 37

    zone_letter = UTMLetterDesignator(lat)

    lon_origin = (zone_number - 1) * 6 - 180 + 3

    lat_rad = math.radians(lat)
    lon_rad = math.radians(lon)
    lon_origin_rad = math.radians(lon_origin)

    ecc_prime_squared = UTM_E2 / (1 - UTM_E2)

    n = WGS84_A / math.sqrt(1 - UTM_E2 * math.sin(lat_rad) ** 2)
    t = math.tan(lat_rad) ** 2
    c = ecc_prime_squared * math.cos(lat_rad) ** 2
    a = (lon_rad - lon_origin_rad) * math.cos(lat_rad)

    m = WGS84_A * (
            (1 - UTM_E2 / 4 - 3 * UTM_E2 ** 2 / 64 - 5 * UTM_E2 ** 3 / 256) * lat_rad
            - (3 * UTM_E2 / 8 + 3 * UTM_E2 ** 2 / 32 + 45 * UTM_E2 ** 3 / 1024) * math.sin(2 * lat_rad)
            + (15 * UTM_E2 ** 2 / 256 + 45 * UTM_E2 ** 3 / 1024) * math.sin(4 * lat_rad)
            - (35 * UTM_E2 ** 3 / 3072) * math.sin(6 * lat_rad)
    )

    easting = (
            UTM_FE
            + UTM_K0
            * n
            * (
                    a
                    + (1 - t + c) * a ** 3 / 6
                    + (5 - 18 * t + t ** 2 + 72 * c - 58 * ecc_prime_squared) * a ** 5 / 120
            )
    )
    if lat < 0:
        northing = UTM_FN_S + UTM_K0 * (m + n * math.tan(lat_rad) * (
                    a ** 2 / 2 + (5 - t + 9 * c + 4 * c ** 2) * a ** 4 / 24 + (
                        61 - 58 * t + t ** 2 + 600 * c - 330 * ecc_prime_squared) * a ** 6 / 720))
    else:
        northing = UTM_FN_N + UTM_K0 * (m + n * math.tan(lat_rad) * (
                    a ** 2 / 2 + (5 - t + 9 * c + 4 * c ** 2) * a ** 4 / 24 + (
                        61 - 58 * t + t ** 2 + 600 * c - 330 * ecc_prime_squared) * a ** 6 / 720))

    return easting, northing, zone_number, zone_letter


def UTMLetterDesignator(lat):
    if 84 >= lat >= 72:
        return 'X'
    elif 72 > lat >= 64:
        return 'W'
    elif 64 > lat >= 56:
        return 'V'
    elif 56 > lat >= 48:
        return 'U'
    elif 48 > lat >= 40:
        return 'T'
    elif 40 > lat >= 32:
        return 'S'
    elif 32 > lat >= 24:
        return 'R'
    elif 24 > lat >= 16:
        return 'Q'
    elif 16 > lat >= 8:
        return 'P'
    elif 8 > lat >= 0:
        return 'N'
    elif 0 > lat >= -8:
        return 'M'
    elif -8 > lat >= -16:
        return 'L'
    elif -16 > lat >= -24:
        return 'K'
    elif -24 > lat >= -32:
        return 'J'
    elif -32 > lat >= -40:
        return 'H'
    elif -40 > lat >= -48:
        return 'G'
    elif -48 > lat >= -56:
        return 'F'
    elif -56 > lat >= -64:
        return 'E'
    elif -64 > lat >= -72:
        return 'D'
    elif -72 > lat >= -80:
        return 'C'
    else:
        return 'Z'  # 错误标志，纬度超出UTM 84N到80S的范围
